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The  formulation  of  a  fully  compatible  general  quadrilateral 
plate  bending  element  is  described.  The  element  is  assembled  from 
four  partially  constrained  linear  curvature  compatible  triangles, 
arranged  so  that  no  mid-side  nodes  occur  on  the  external  edges  of 
the  quadrilateral;  thus,  the  resulting  element  has  only  12  degrees 
of  freedom.  Also  described  is  a  simple  shear  distortion  mechanism 
which  may  be  incorporated  into  the  element  without  changing  its 
basic  structure.  Results  are  presented  for  static  analyses  with  and 
without  shear  distortion,  and  for  plate  vibration  and  plate  buckling 
studies,  all  performed  with  this  quadrilateral  element.  It  is  con¬ 
cluded  that  this  is  the  most  efficient  general  bending  element  yet 
devised. 
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SECTION  I 
INTRODUCTION 

HISTORICAL  BACKGROUND 

The  finite  element  method  can  be  considered  as  the  most  powerful  and  versatile  dis¬ 
cretization  technique  presently  available  for  the  numerical  solution  of  complex  structural 
problems  using  digital  computers.  The  method  was  developed  originally  as  an  application 
of  standard  structural  analysis  procedures  to  a  physically  discretized  approximation  of  the 
actual  system;  the  concept  has  been  extensively  described  elsewhere  (References  1,  2,  3, 
and  4)  and  will  not  be  detailed  here.  During  the  past  few  years,  study  of  the  mathematical 
foundations  of  the  method  (References  5,  6,  7,  and  8)  as  well  as  its  application  to  a  wider 
class  of  field  problems  (Reference  4)  has  greatly  clarified  the  basic  requirements  for  its 
effective  formulation. 


The  application  of  the  finite  element  method  to  plate  bending  problems  dates  from  the 
late  1950’s.  The  first  successful  results  were  published  in  1960  (References  9,  10,  and  11), 
but  were  obtained  using  rectangular  elements,  none  of  which  satisfied  the  requirements 
listed  in  the  following  section.  The  construction  of  adequate  displacement  expansions  for 
the  more  versatile  triangular  and  quadrilateral  shapes  was  not  achieved  until  1965,  most 
of  these  results  being  presented  at  the  1st  Air  Force  Institute  of  Technology  (AFIT)  Con¬ 
ference  (References  12,  13,  and  14).  By  this  time,  basic  guides  for  the  selection  of  suitable 
element  deformation  patterns  had  been  set  forth  empirically  (Reference  15)  and  later 
rigorously  proved  (References  7,  and  8).  Most  of  the  early  developments  in  plate  bending 
analysis  by  finite  elements  are  summarized  in  Reference  12.  During  the  three  years  since 
the  1st  AFIT  Conference,  a  considerable  number  of  publications  has  been  devoted  to  general 
plate  bending  elements.  Several  displacement-assumed  models  with  various  degrees  of 
refinement  have  been  proposed  (References  16,  17,  18,19,  20,  21,  and  22)  as  well  as  elements 
for  analyses  based  on  equilibrium  (References  14,  23,  and  24),  mixed  (References  25,  26, 
27,  28,  and  29)  and  “hybrid”  (References  30,  and  31)  variational  principles. 


In  view  of  this  long  list  of  alternatives,  it  is  important  to  consider  which  alternative 
provides  the  best  balance  in  practical  usage,  taking  into  consideration  factors  such  as 
simplicity  of  formulation,  versatility  of  application,  reliability,  computational  effort,  and 
accuracy.  The  purpose  of  this  paper  is  to  describe  a  general  quadrilateral  element,  formed 
as  an  assemblage  of  triangular  elements,  which  is  believed  to  be  one  of  the  most  efficient 
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mesh  units  for  both  plate  and  thin  shell  applications.  The  quadrilateral  element  is  designated 
Q-19  to  identify  it  as  a  quadrilateral  with  19  basic  degrees  of  freedom.  (It  is  reduced  to 
12  degrees  of  freedom  before  incorporating  into  an  element  assemblage.)  The  triangular 
elements  of  which  it  is  formed  are  designated  LCCT-11,  meaning  Linear  Curvature  Com¬ 
patible  Triangle  with  11  degrees  of  freedom.  The  9-degree  of  freedom  version  of  this 
element  (LCCT-9)  is  the  same  as  that  designated  HCT  in  Reference  12. 

The  element  stiffness  derivation  is  outlined  completely,  after  the  basic  requirements  and 
limitations  imposed  on  the  assumed  displacement  expansions  are  summarized.  Then  results 
obtained  in  various  static,  dynamic,  and  buckling  analyses  are  presented  to  demonstrate  the 
effectiveness  of  the  element. 

DISPLACEMENT  ASSUMPTION  REQUIREMENTS 

The  transverse  displacement  w(x,y)  is  the  sole  primary  variable  required  in  the  formula¬ 
tion  of  the  total  potential  energy  of  a  plate  element.  The  set  of  assumed  displacement  patterns 
in  each  element  should  belong  to  the  admissible  class  of  functions  satisfying  the  following 
requirements: 

(1)  Compatibility:  (a)  The  assumed  w(x,y)  mustbe  continuous  and  have  continuous  first 
derivatives  inside  each  element;  (b)  w(x,y)  and  its  normal  derivatives  dw/dn  (the  normal  slope) 
must  be  uniquely  specified  along  any  element  interface  S  (where  n  is  normal  to  S)  by  nodal 
displacement  selected  on  S. 

(2)  Completeness:  All  rigid  body  displacement  states  and  uniform  strain  (constant 

2  2 

curvature)  states  must  be  included  in  the  expansion.  In  other  words,  the  six  terms  l,x,y,x  ,xy,y 
(or  their  equivalent  in  other  coordinate  systems)  must  be  included  in  the  set  of  element 
displacement  modes. 

If  these  requirements  are  met,  the  finite  element  analysis  is  a  special  form  of  the 
classical  Rayleigh-Ritz  procedure  (Reference  5)  in  which  the  nodal  displacements  are  taken 
as  generalized  coordinates.  For  a  stable  elastic  material,  it  follows  from  the  positive 
definite  character  of  the  total  potential  energy  functional  that  a  sequence  of  finite  element 
solutions  obtained  by  a  mesh  subdivision  process  provides  a  minimizing  sequence  for  the 
strain  energy.  Using  Sobolev’s  inclusion  theorem  (Reference  32),  it  can  be  shown  that  the 
transverse  displacements  converge  uniformly  whereas  slopes,  curvatures  and  bending 
moments  converge  in  the  mean.  (If  the  exact  solution  is  sufficiently  smooth,  rotations  and 
curvatures  also  converge  uniformly.) 
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The  compatibility  requirement  is  not  actually  necessary  for  strain  energy  convergence. 
“Non-conforming”  elements  which  violate  slope  continuity  between  corners  have  been 
used  extensively  and  often  with  good  results  (References  5,  10,  11,  12,  13,  and  24).  If  any 
group  of  non-conforming  elements  can  represent  rigid  body  and  constant  curvature  states 
exactly,  the  finite  element  solutions  will  converge  to  the  true  value  of  strain  energy  (Ref¬ 
erence  7) ,  though  not  necessarily  in  monotonic  form.  Although  their  derivation  is  simpler, 
the  use  of  non-conforming  elements  in  general  purpose  programs  has  certain  disadvantages: 
(a)  energy  convergence  depends  on  the  mesh  subdivision  pattern  (Reference  13),  (b)  curvatures 
and  bending  moments  may  not  converge  even  if  the  strain  energy  does,  and  (c)  no  error 
control  is  available. 

CONSTRUCTION  OF  DISPLACEMENT  EXPANSIONS 


In  considering  the  selection  of  transverse  displacement  patterns  for  general  polygonal 
flat  plate  elements,  it  is  convenient  to  distinguish  between  two  classes: 

2 

1.  Class  C  :  the  assumed  w(x,y)  has  continuous  second  derivatives  (curvatures)  at 
the  element  corners  (inside  the  element).  Example:  single  polynomial  expansions. 

2.  Class  C1:  w(x,y)  may  have  discontinuous  corner  curvatures. 

2 

The  generation  of  plate  elements  using  C  expansions  is  severely  restricted  by  the 
following  limitation  principle  (which  is  proved  in  Reference  18): 

2 

To  construct  a  complete  and  compatible  C  plate  bending  expansion,  a  minimum  of 

6  degrees  of  freedom  (w>wx,wy'wxx>wxy  an<^  w^^)  are  required  at  each  non-right  angled 
corner,  and 


4  degrees  of  freedom  (w,w  ,w  ,w  , 

x  y  xy 

required  at  each  right-angled  corner. 


where  x-y  are  taken  along  the  adjacent  sides)  are 


If  less  degrees  of  freedom  are  selected  and  full  compatibility  is  enforced,  completeness  is 

lost  and  the  finite  element  analysis  does  not  converge.  It  follows  that  at  least  18  degrees  of 

freedom  are  required  for  a  general  triangle,  24  for  an  arbitrary  quadrilateral  and  16  for  a 

rectangle.  When  polynomial  expansions  are  used,  at  least  a  quintic  (which  provides  21  gen- 

2 

eralized  coordinates)  is  necessary  for  a  C  expansion  over  a  triangle  (References  18  and  19.) 
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The  use  of  C1  expansions  provides  more  flexibility  in  the  selection  of  displacement 
functions,  and  is  essential  for  generating  compatible,  complete  elements  with  only  three 
degrees  of  freedom  (w.w^w^)  per  comer.  Two  construction  methods  have  been  used  with 
success: 

(1)  C1  spline  fit  of  several  polynomial  subexpansions  assumed  over  triangular  sub- 
regions  (References  12,  14  and  20).  (The  element  described  in  this  paper  fits  into  this 
category.) 

(2)  Correction  of  single  polynomial  expansions  with  rational  functions  exhibiting 
curvature  singularities  at  the  corners  (References  13  and  20). 

Another  limitation  principle  concerning  the  construction  of  curved  side  plate  elements 
is  given  in  Reference  8. 


SECTION  n 

DERIVATION  OF  THE  ELEMENT  STIFFNESS  MATRICES 

STIFFNESS  OF  THE  LCCT-12  ELEMENT 
Triangle  Geometry 

The  geometry  of  an  arbitrary  triangular  element  can  be  expressed  in  a  Cartesian 
coordinate  system  by  its  nodal  coordinates  or  its  projected  dimensions,  as  shown  in  Figure  la, 
or  alternatively  by  its  intrinsic  dimensions  as  defined  in  Figure  lb.  If  j  and  k  denote  the 
first  and  second  cyclic  permutations  of  i  =  1,2,3  (i.e.,  j  =  2,3,1  and  k  =  3,1,2),  the  projected 
dimensions  may  be  expressed 

«i  ■  *k-*j  :  bi  *  »r»i.  i" 


403 


AFFDL-TR-68-150 


Also  the  intrinsic  dimensions  may  be  defined  in  terms  of  the  projected  dimensions,  one  of 
the  more  important  relationships  being 

dj  =-(ajOk+bi  bk)/Lj  (2) 


The  analysis  of  the  stiffness  properties  of  a  triangular  element  is  greatly  simplified  by 
the  use  of  triangular  (natural)  coordinates.  The  triangular  coordinates  £2>  £3  of  any 
point  “P”  in  the  triangle  may  be  defined  either  as  the  ratios  of  the  areas  A^  of  the  sub- 
triangles  subtended  by  that  point  to  the  total  area  A  of  the  triangle,  or  as  the  ratios  of  the 
normal  distances  n.  to  the  heights  H.;  i.e. 


as  shown  in  Figure  2.  It  should  be  noted  that  the  triangular  coordinates  are  related  by  the 
constraining  condition  ^  +  ^2  +  ^3  =  *" 


The  relationships  between  Cartesian  and  triangular  coordinates  may  be  expressed  as 
follows 


or  by  inversion: 


1 

1  1  1 

c, 

X 

= 

X 

CM 

X 

X 

y 

-  yi  y2  y 3 . 

_c3. 

c2 

1 

=  2  A 

_2A23 

2  A  31 

b|  °l 

b  2  a2 

1 

X 

c. 

ro 

> 

b  3  °3 

>1 

- _ _ 1 

(4) 


(5) 


in  which  2A.  ,  =  x.y..  -  x^y..  From  Equation  5,  the  following  important  differential  relationships 
between  the  coordinate  systems  may  be  noted: 


^Cj  _  _!j__  a£j  .  gi 

dx  '  2A  dy  2  A 


(6) 
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Subelement  Geometry 

It  was  statedinthelntroductionthatitis  impossible  to  derive  a  fully  compatible  triangular 
plate  bending  element  using  only  a  single  cubic  displacement  expansion.  To  obtain  a  compatible 
system,  the  basic  element  has  been  divided  into  three  subelements,  as  shown  in  Figure  3a  in 
which  the  internal  point  “O'*  is  the  centroid  of  the  element  area,  and  the  subtriangles  are 
numbered  to  correspond  with  the  opposite  corner  number.  The  geometric  relationships 
discussed  above  apply  to  each  of  the  subtriangles  if  its  three  comers  are  renumbered  1-2-3, 
with  3  being  taken  as  the  internal  point.  The  subtriangle  number  is  then  identified  in  the 
algebraic  expressions  by  means  of  a  superscript.  The  renumbering  scheme  for  Subelement 
1  is  shown  in  Figure  3b. 

Displacement  Interpolation  Functions 

The  nodal  displacement  degrees  of  freedom  which  are  to  be  considered  in  the  stiffness 
matrix  of  the  complete  element  also  are  shown  in  Figure  3a,  These  include  the  transverse 
displacements  of  each  corner,  w.,  the  rotations  at  each  corner  about  the  x  and  y  axes, 
#x.  and  9 as  well  as  the  rotation  at  the  three  mid-side  nodes  about  an  axis  parallel  to 
the  side,  0^,  0_  and  0g.  In  order  to  develop  displacement  patterns  for  the  complete  triangle, 
displacement  interpolation  functions  were  assumed  independently  for  each  subtriangle. 

The  displacement  interpolation  functions  for  each  subelement  express  the  relationship 
between  the  displacements  w^  within  the  element  and  the  ten  displacement  components  of 
its  nodal  points  r  ^ .  as  follows : 

«<*>=  17) 

As  may  be  noted  in  Figure  3b,  the  nodal  displacement  vector  for  subelement  1  is: 

r  [  w2  w3  ^X3  &yz  wo  ^xo  ^yo  ^5  ] 
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Figure  2.  Triangular  Coordinate  System 


(a)  Element  Assembly  (b)  Subelement  Nodal  System 

Figure  3.  Assembly  of  the  LCCT-12  Plate  Bending  Element 
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The  set  of  10  cubic  interpolation  functions  for  Subelement  1  may  be  expressed  in  triangular 
coordinates  as  follows: 


C2  ( 3-  2  £  )+  6^ll)  £  £  £ 

I  H  *3  b2  =3 


Vb3  St-bz5j)f(bs  Mrb,  >  C,C,C 


I  •’2 33 


r2,  (l)r  (l)  r  ..  ,  (l)  (l)  (0,  t  r 

C,  °  3  C2  "a2  C3)+(°3  -0,  )  C,C2C, 

Cj  (3-  2  £2)+  6X(^  £,£.£. 


£2  (bt‘)£  -b(,)£  )+  £,  t  £, 

5  2  ls3  31  2  3  3  123 

2  (l).  (•)„  (')  (,)ll)  „  „  „ 

^2^°l  C3  °3  C  |  ^  2  "  °3  ^3*  C|  C2  C3 

C|  <3-2C3  ) 


.2 
•3 

2  (I)  „  (I)  „ 

C3  (bVCr»V  fit  ) 


r2  .  (i)  r  I'L  > 

^3  °  2  C,-a,C2) 
C,C3 


(9) 


where  the  subscripts  correspond  to  the  renumbered  nodes  of  the  subelement.  With  this 
convention,  the  interpolation  functions  for  subelements  2  and  3  are  the  same  as  Equation  9, 
except  with  appropriately  changed  superscripts.  It  should  be  noted,  however,  that  the  nodal 
displacements  in  Equation  8  are  identified  by  node  numbers  defined  for  the  complete  element 
assembly. 


Now  if  the  vector  r  of  all  nodal  displacements  of  the  complete  element  assembly  are 
written  in  the  sequence 


-T 

r 


[  WI  ®xt  ®yi  w2  ^X2  ^y  2  w3®x3^y3  ^4  ^siwo^xo^yo 


■  W  I  '0  ] 


UO) 


the  displacements  in  subelement  1  can  be  expressed 
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where  <£''*  is  similar  to  Equation  9,  but  expanded  with  5  zeros  to  account  for  the  nodal 
displacements  not  associated  with  Element  1,  and  with  appropriate  rearrangement  of  terms. 
The  submatrices  *  and  ^  represent  the  interpolation  functions  for  the  external  and 
internal  nodal  displacements  respectively. 

Expressing  the  displacements  in  the  other  subelements  similarly,  the  complete  system 
of  displacements  can  be  written: 


Establishing  Internal  Compatibility 

Equation  12  is  an  expression  of  the  cubic  displacement  patterns  developed  in  the  three 
subelements.  Because  of  the  common  displacements  imposed  at  the  nodes,  the  transverse 
displacements  of  two  adjacent  elements  are  identical  along  their  juncture  line.  However, 
their  normal  slopes  differ  between  the  nodes;  hence.  Equation  12  does  not  represent  an 
internally  compatible  displacement  field.  To  establish  slope  compatibility  along  the  internal 
edges  of  the  subelements,  additional  nodes  7,  8,  and  9  were  located  at  the  mid -points  of 
these  edges,  as  shown  in  Figure  3a.  The  normal  slope  was  computed  at  each  of  these  nodes 
in  each  subelement,  for  example 
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Equation.  14  may  be  expressed  symbolically  as  follows: 

tB  !B» ][>;-]=  0 

Now  the  values  of  tq  which  will  satisfy  these  compatibility  conditions  may  be  computed 
from  Equation  15;  i.e. 

ro  s  —  Bo  B  r  =  L  r  (16) 

Finally,  introducing  the  compatibility  constraint  of  Equation  16  into  Equation  12,  the  fully 
compatible  displacement  field  in  the  three  subelements  becomes: 


[  *(1)1 
M) 

\ 

= 

+ 

<*(z) 

M) 

L 

r  = 

v/3) 

lO 

"V 

_ i 

_~0  . 

Although  it  is  straightforward  in  concept,  the  derivation  of  the  compatible  interpolation 
functions  of  Equation  17  involves  long  algebraic  manipulations  and  is  a  tedious  operation. 
Explicit  expressions  for  these  functions  are  presented  in  the  Appendix  for  the  convenience 
of  the  reader. 


Internal  Curvature  Field 


In  order  to  define  the  stiffness  matrix  corresponding  with  the  derived  displacement 
field,  it  is  necessary  to  establish  the  curvatures  within  the  subelements.  The  most  convenient 
means  of  defining  the  curvature  distribution  is  by  means  of  appropriate  interpolation  functions; 
in  the  present  case  where  the  assumed  displacement  functions  are  cubic,  it  is  important  to 
note  that  the  curvature  must  vary  linearly  within  each  subelement.  Thus,  the  curvature 
field  can  be  expressed  as  the  product  of  linear  interpolation  functions  multiplied  by  curvature 
values  defined  at  the  corner  nodes. 


Within  any  subelement  “i,”  the  curvature 
the  displacement  field,  thus 


can  be  obtained  by  differentiation  of 


(i  ) 
wxx 


w 

2W 


(i) 
yy 
(i ) 
xy 


T(i)r 


(18) 
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where 


T(i  = 


ax2 

a2*(i> 

ay2 

dz4 


*£(i) 


L  ax  ay  J 


(19) 


Now  the  nodal  values  of  the  curvatures  in  subelement  “i”  may  be  determined  by  evaluating 


T(l>  at  corner  points;  thus,  if  the  nodal  curvatures  are  designated  and  the  nodal 

(i>  _  (i)  n 


(i) 


values  of  T 


are  T 


the  relationship  may  be  expressed 


(i ) 


<n  r 


(20) 


in  which,  for  example, 


'(i)T  =  f  w(i)  w(i)  wU)  W(i)  w(i)  w(i)  2w(i>  2wli)  2*(i)  1 

-n  L  xxz  wxx3  wxxo  "y y2  wyy3  wyyo  ^wxyz  ^wxy3  ^wxyo  J 


(21) 


The  linear  curvature  variation  within  the  subelement  can  now  be  expressed  by  means  of  the 
linear  interpolation  functions  x  as  follows: 


,( i)  . 


*0  0 


X 

o  *x 
0  0 


0 

f. 


v{i) 

*  n 


-  *x  X  n 


(i) 


(22) 


where  the  linear  interpolation  functions  are  merely  the  three  triangular  coordinates: 

*x  =  [  C,  C2  C3  ]  (23) 

It  is  of  interest  to  note  that  the  curvatures  in  the  three  subelements  are  identical  at  the 
nodal  point  “O”;  thus,  these  quantities  need  be  evaluated  only  once. 
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Element  Stiffness 


The  strain  energy  due  to  bending  of  the  subelement  may  be  obtained  by  integrating  the 
product  of  the  moments  and  the  curvatures  over  the  area  of  the  element,  i.e. 

U(i)  fA  m(i)T  X(i>dA  (24) 

where  the  internal  moments  are  given  by 


=  DX  (i)  (25) 


m 


(  i)  = 


m 


m 


m 


(i) 

xx 

(i) 

yy 

U) 

xy 


The  matrix  D  in  Equation  25  represents  the  constitutive  relationship  for  the  element 
material.  For  a  plate  with  constant  material  properties  through  the  thickness  “h,”  this 
relationship  may  be  expressed 


D 


(Symm.) 


(26) 


where  the  coefficients  C. ^  represent  the  elastic  properties.  Substituting  Equations  20,  22  and 
25  into  Equation  24  leads  to 


dA  T 


(i) 


r 


or 


K^r 


(27) 


Where 


(28) 
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is  the  stiffness  matrix  contribution  of  Subelement  “i”  expressed  in  terms  of  the  nodal 
Displacements  r  ,  and 


(i) 

G 


r(i) 
G  12 
,(i) 
>22 


(Symm.) 


(i 

'23 

(i 

3 


(29) 


If  the  material  properties  and  thickness  are  uniform  over  the  area  of  the  subelement,  the 
individual  terms  of  the  Matrix  G  ^  become 


where 


R 


(i) 


is  given  by 


(i) 


£  f  ^  =  V 


2  I 
I  2 
I  I 


I 

1 

2 


(30) 


If  material  properties  or  thickness  of  the  element  vary  over  its  area,  it  is  recommended  to 
evaluate  the  intergral  of  Equation  29  numerically. 


It  will  be  noted  that  the  stiffness  matrix  of  the  complete  triangle  is  obtained  by  merely 
adding  the  contributions  of  the  three  subelements  because  they  are  all  expressed  in  terms 
of  the  same  set  of  nodal  coordinates.  Thus 


K 


K('  1  + 


.(2) 


,(3 


(31) 


ELEMENT  STIFFNESS  INCLUDING  SHEAR  DISTORTION 


Assumed  Displacements 


The  basic  assumption  of  the  element  stiffness  analysis  presented  in  the  preceding 
paragraphs  is  the  Kirchhoff  hypothesis,  which  maybe  expressed  mathematically  as  follows; 


w  (  x,  y,  z  )  =  w(x,y) 


u 

V 


z 


z 


d  w 
dx 
d  w 
dy 


> 


/ 


(32) 
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This  assumption  imposes  the  condition  that  normals  to  the  undeformed  plate  mid-surface 
remain  undistorted  and  normal  to  the  deformed  mid-surface.  The  fact  that  they  remain 
normal  during  deformation  effectively  eliminates  shear  distortion  from  the  assumed  behavior, 
and  thus  leads  to  results  which  are  applicable  only  to  relatively  thin  plates.  The  theory 
can  be  extended  to  account  for  shear  distortion  in  an  approximate  way,  however,  by  adding 
a  simple  shear  distortion  mechanism  to  the  Kirchhoff  deformation  hypothesis. 


The  shear  distortion  mechanism  assumed  in  this  study  can  be  explained  conveniently  by 
reference  to  Figure  4a,  which  represents  a  cross-section  view  of  a  plate  element.  The 
rotation  of  the  cross-section  is  shown  to  depend  on  the  rotation  of  the  mid-surface  d7plus 
an  additional  shear  distortion  /3  which  is  assumed  to  be  a  simple  straight  line  rotation 
(uniform  shear  strain  through  the  thickness).  The  total  rotation  <£  K  thus  is  given  by 


d  w 
dy 


(33) 


A  similar  assumption  is  made  along  the  other  axis,  thus 

dw 
d  x 


<#>,=  — -0, 


By  -0, 


(34) 


To  define  this  distribution  of  cross-section  rotations  through  the  element  field,  it  is 
now  assumed  that  the  transverse  displacements  w(x,y)  are  given  by  the  compatible  inter¬ 
polation  functions  defined  by  Equation  17,  for  each  subelement: 


»(i)  B 


(35) 


(a)  Definition  of  Mean  Shear  Distortion 


(b)  Continuity  Requirement 


Figure  4.  Shear  Distortion  Geometry 
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In  addition,  it  is  assumed  that  the  shear  distortions  vary  linearly  over  the  entire  triangle: 


or  symbolically 


c, 


C2 


£3  0 

0  c, 


c2 


£3 


0*. 

0X2 

0X3 

0yi 

0y2 

0y  3 


*  '  */3  '» 


(36) 


(37) 


Now  it  is  convenient  to  substitute  total  rotation  of  the  section  <£  for  the  slope  Q  in  the 
nodal  degrees  of  freedom  of  Equation  35.  Thus,  using  relationships  derived  from  Equation 
33  and  34,  the  following  expressions  may  be  defined  for  each  comer  nodal  point: 


re,i' 

£xi 

-0yi 

CD 

*■< 

1 _ 

_4>  yi 

_+0xi 

while  for  a  typical  mid-side  node,  such  as  node  4: 

04  =  <£4  -  <0X1  +0X2>T  "  {0yi+0y2)-f  (39> 

where  Cg  =  ag/Lg  and  Sg  =  -bg/Lg.  Similar  expressions  may  be  obtained  for  nodes  5  and  6 
by  cyclic  permutation  of  indices. 


Using  these  expressions,  the  interpolation  functions  for  displacement  in  each  subelement 
Equation  35  may  be  rewritten  as  follows: 


where 


[  ",  <f>yt  w  2  0X2^2  W  <£5tf>6  ] 

displacement  vector  in  Eq.36 
interpolation  functions  of  Eq,  17 


(40) 
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and 


~s 


'^0y. 

-  c3  +  fye  c2>  /z 

^0y2 

'  l*6'  * +B,  C*)/2 

“*0y» 

°2  °i)/2 

-*0k, 

-%,s>+*8'  s‘l/2 

^0X2 

S'  +  ^4  *»’  /2 

^0X3 

~(c^0  6  S2+<£05S|)/2 

(41) 


in  which  the  subscripts  of  the  interpolation  functions  4>  identify  their  nodal  displacement 
components.  Equations  36,  37  and  40  define  completely  the  deformation  of  the  element, 
expressing  the  transverse  displacements  Equation  40  in  terms  of  18  nodal  displacements, 
and  the  shearing  distortions  Equation  36  in  terms  only  of  the  6  nodal  shear  distortions. 

Internal  Curvature  Field 

Where  shearing  distortions  are  included  in  the  displacement  field,  the  total  curvatures 
X*  which  define  the  normal  strain  distribution  depend  on  both  the  transverse  displacements 
and  the  shear  distortions  as  follows: 


dz  w 

d  x 

dx2 

dx 

dy 

r 

d2w 

<3y2 

-+ 

a/3, 

dy 

dcf>  d(/>x 

*  j 

^2 

0  w 

o  . 

d/3x  d/3 

l  * 

dy  dx 

dx  dy 

dy  dx 

(42) 


The  first  curvature  term  must  be  defined  separately  for  each  subelement,  and  from  Equation 
40  is  given  by: 


(43) 
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where  T  is  given  by  Equation  18  and  T  ^  has  an  equivalent  definition  in  terms  of 
<£(  i )  — 

<PS  •  However,  the  second  curvature  term  in  Equation  42  may  be  expressed  for  the 
entire  element,  i.e.: 


where 


(44) 


(45) 


applying  the  derivative  definitions  of  Equation  6  to  Equation  36.  Thus  substituting  Equations  43 
and  44  into  Equation  42  leads  to 


where 


.(i) 


_rJ?_ 

L  !  s  J 

rs 

L  J 

(46) 


(47) 


Now  because  these  combined  curvatures  still  vary  linearly  within  the  subelements,  their 
distribution  can  be  expressed  in  terms  of  nodal  curvature  values  as  before;  i.e.  by  analogy 
with  Equations  20  and  22: 


and 


Element  Stiffness 


(48) 


(49) 


The  strain  energy  due  to  bending  deformations  may  be  expressed  in  terms  of  the  bending 
stiffness  matrix  by  analogy  with  Equation  31: 


(  i) 

B 


"ffirr 

sn 


(i) 

sn 


] 


(50) 
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in  which  q*'*  is  exactly  as  in  Equation  30.  An  additional  strain  energy  contribution  Ug 
results  from  the  shear  distortion,  however,  which  is  given  as  follows: 


U 


s 


in  which  the  shearing  forces  Q  are 


(51) 


(52) 


In  this  expression,  Ds  represents  the  shear  distortion  constitutive  relationship.  For  a 
plate  with  uniform  properties  through  the  thickness,  and  assuming  that  z  is  a  principal 
elastic  direction,  it  may  be  expressed: 


«>s 


h 


C44  C45  "I 
C43  ^59 


in  which  C..  are  material  constants, 
i] 

Introducing  Equations  53,  52  and  37  into  Equation  51  leads  to 


(53) 


’  "2  'll* 8Ds*/3  dAr» 


or 


K  r 

ss  s 


(54) 


where  K“  is  the  shear  stiffness  matrix.  In  the  case  where  material  properties  and  thickness 
ss 

are  uniform  over  the  element,  the  shear  stiffness  becomes 


K 


1 1 

ss 


r  C44  R  C45  R 
c  45  R  C55  R 


(55) 


in  which  R  is  given  by  Equation  30.  For  nonuniform  thickness  or  material  properties, 
the  integral  should  be  evaluated  numerically. 


The  total  stiffness  is  obtained  finally  by  superposition  of  the  bending  and  shear 
contributions, 


k_bb  I  _  R  BS_ _ 

KSB  j(Kss  +  KSS  ] 


(56) 
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in  which  it  should  be  noted  that  the  bending  terms  must  be  computed  separately  for  each 
subelement  and  added,  whereas  the  shear  stiffness  term  K"  may  be  computed  for  the 

O  V 

entire  element. 


Condensation  of  Shear  Distortion  Degrees  of  Freedom 


The  element  stiffness  matrix  of  Equation  56  represents  the  element  force-displacement 
relationship,  which  may  be  expressed  as  follows: 


PB 

K  BB 

KBS 

“ 

rs 

Si 

1 - 

“0 

i  10 

K  SB 

KSS 

t 

CO 

1 _ 

in  which  Kcc=  K'  +  K .  A  total  of  18  degrees  of  freedom  are  included  in  this  re- 
lationship,  6  shear  distortion  components  rs  in  addition  to  the  12  basic  translation  and 
rotation  displacements  of  the  nodes.  At  this  point,  however,  it  should  be  noted  that  the 
shear  distortion  angle  need  not  be  continuous  between  adjacent  elements.  As  shown  in 
Figure  4b,  the  total  rotation  <f>  must  be  the  same  for  two  adjacent  elements  in  order  to 
satisfy  compatibility  conditions,  but  J3  and  9  need  not  be  the  same.  Consequently,  the  shear 
distortion  degrees  of  freedom  are  not  needed  to  achieve  compatibility  and  they  may  be 
eliminated  from  the  element  by  static  condensation. 


As  a  result  of  the  condensation  process,  the  modified  element  stiffness  relationship 
becomes: 


P  z  K  r 


(58) 


in  which 


K 

P 


Kss  ksb 

K-1  P 

*ss  s 


(59) 


It  is  of  interest  to  note  that  the  degrees  of  freedom  of  this  reduced  element  stiffness  matrix 
are  exactly  equivalent  to  those  of  the  original  LCCT-12  element;  thus,  an  element  including 
this  shear  distortion  capability  can  easily  be  incorporated  into  existing  plate  bending  analysis 
programs  based  on  the  original  element. 
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GENERAL  QUADRILATERAL  ELEMENT  Q-19 

Although  the  LCCT-12  element  employs  an  optimum  compatible  cubic  displacement 
field  and,  therefore,  will  yield  the  best  possible  results  for  a  given  triangular  element 
mesh  involving  compatible  cubic  displacements,  its  mid-point  nodes  are  a  somewhat  un¬ 
desirable  feature.  They  tend  to  complicate  mesh  generation  procedures,  increase  the  band 
width  of  the  assembled  equation  systems,  and  require  special  identification  in  the  development 
of  computer  programs.  To  overcome  these  disadvantages,  while  retaining  most  of  the 
flexibility  of  the  LCCT-12  element,  it  is  convenient  to  develop  a  special  version  of  the  element 
by  constraining  the  normal  slope  to  vary  linearly  along  one  side. 


Consider,  for  example.  Subelement  3  of  the  element  shown  in  Figure  3a.  The  mid-side 
node  4  of  this  subelement  can  be  eliminated  by  introducing  the  condition  that  its  value  is 
the  average  of  the  corresponding  slopes  at  nodes  1  and  2, 

0,  =  0x,  c3  -  0,«  s3  =  <  ’T“  1601 

where  c3  and  s3  are  the  same  as  in  Equation  39.  Using  this  condition,  the  displacement 
interpolation  expressions  in  Equation  17  can  be  reduced  to  only  11  components,  and  the 
stiffness  matrix  thus  reduced  accordingly. 


The  resulting  partially  constrained  element  is  designated  LCCT-11.  (Similar  con¬ 
straints  must  be  applied  to  the  other  sides  to  develop  the  LCCT-10  and  LCCT-9  elements.) 
Four  LCCT-11  elements  may  then  be  assembled  into  the  Q-19  quandrilateral  having  no 
mid-side  nodes  on  the  exterior  edges  as  shown  in  Figure  5.  Although  this  element  has 
19  degrees  of  freedom,  the  7  internal  degrees  of  freedom  of  the  assemblage  are  eliminated 
by  a  static  condensation  process  equivalent  to  Equation  56  before  the  quadrilateral  is  as¬ 
sembled  into  the  complete  structure.  Thus,  the  final  quadrilateral  has  only  12  degrees  of 
freedom,  corresponding  to  a  translation  and  two  rotations  at  each  node.  It  is  a  fully  com¬ 
patible  element,  having  linear  variations  of  normal  slopes  along  all  exterior  edges. 

Also  it  should  be  noted  that  triangular  elements  including  shear  distortion  effects  can 
be  assembled  to  form  a  quadrilateral  element  in  a  similar  fashion.  In  this  case,  however, 
it  has  been  found  most  convenient  to  condense  the  shear  distortion  degrees  of  freedom  at  the 
quadrilateral  assemblage  level  rather  than  in  the  individual  triangles. 
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3 


Figure  5.  Quadrilateral  Q-19  Assembled  with  4  LCCT-11  Elements 


SECTION  in 

RESULTS  OF  ANALYSES  WITH  THE  Q-19  ELEMENT 

STATIC  PLATE  BENDING 

Neglecting  Shear 

In  order  to  evaluate  the  accuracy  obtainable  with  the  Q-19  quadrilateral  element,  a 
series  of  convergence  studies  were  made  of  the  central  deflection  developed  in  a  square 
plate.  Both  uniform  pressure  and  concentrated  central  loading  cases  were  considered; 
boundary  conditions  were  both  simply  supported  and  clamped.  Because  of  double  symmetry, 
only  one-quarter  of  the  plate  was  considered  in  the  analyses. 

The  computed  central  deflections  are  presented  in  Figure  6,  the  abscissa  in  the  graphs 
representing  the  number  of  elements  along  one  side  of  the  square  mesh.  Results  for  the 
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(a)  Central  Deflection  w  =  a  q  a  /D 
of  Simply  Supported  Square  Plate 
Under  Uniform  Load  q 


(c)  Central  Deflection  w  =  P  a  /D 
of  Simply  Supported  Square  Plate 
Under  Concentrated  Central  Load 
P 


4 

(b)  Central  Deflection  w  =  a  q  a  /D 
of  Square  Clamped  Plate  Under 
Uniform  Load  q 


(d)  Central  Deflection  w  =  /3  P  a  /D 
of  Clamped  Square  Plate  Under 
Concentrated  Central  Load  P 


MEANING  OF  SYMBOLS 

M,  ACM 

Incompatible  but  complete  rectangular  elements  (see  Ref.  12);  12  external  D.  F. 

P 

Compatible  but  incomplete  rectangular  element  (see  Ref.  12);  12  external  D.  F. 

LCCT-9 

Two  LCCT-9  triangles  (labeled  HCT  in  Ref.  12);  12  external  D.  F. 

LCCT-12 

Two  LCCT-12  triangles;  16  external  and  1  internal  D.  F. 

Q-19 

Four  LCCT-11  triangles;  12  external  and  7  internal  D.  F. 

DVS 

Compatible  quadrilateral  from  Refs  14  and  23;  16  external  D.  F. 

l 

1 

1 

1 

1 

Indicates  a  rectangular  element  constructed  with  a  single  C2  expansion. 

Denotes  a  triangular  or  quadrilateral  element  constructed  with  a  C ^  expansion. 

Figure  6.  Convergence  Study  of  Various  Plate  Bending  Elements 
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elements  described  in  this  paper  are  labeled  Q-19,  LCCT-12,  and  LCCT-9.  Also  shown  are 
results  for  rectangular  elements  presented  in  Reference  12  (note  that  LCCT-9  was  labeled 
HCT  there),  as  well  as  a  new  quadrilateral  D VS  from  Reference  23,  The  superiority  of 
Q-19  over  any  other  quadrilateral  element  with  12  degrees  of  freedom  is  clearly  evident  —  the 
improvement  over  LCCT-9  is  quite  marked.  It  may  be  noted  that  DVS  gives  slightly  better 
results  than  Q-19,  but  it  has  16  degrees  of  freedom,  including  mid-side  nodes  which  greatly 
hamper  its  computational  efficiency. 


A  second  static  analysis  example  is  presented  in  Figure  7,  the  rectangular  plate  being 
loaded  uniformly,  simply  supported  at  two  sides,  clamped  and  free  at  the  other  two  edges 
(problem  taken  from  Reference  25).  Results  of  a4  x  4  mesh  of  LCCT-12  elements  (2  triangles 
to  each  square)  are  shown.  Displacements  are  essentially  exact,  while  the  moments  are  very 
close  to  the  exact  results.  Computed  moments  from  an  8  x  8  mesh  are  indistinguishable 
from  the  exact  curves  except  at  the  edges.  It  is  probable  that  the  Q-19  element  would  have 
given  still  better  results  and  would  have  required  less  computer  effort.  (However,  the 
8x8  LCCT-12  analysis  required  only  48  seconds  on  an  IBM  7094  computer.) 


Including  Shear  Distortion 


Results  of  the  first  test  of  the  shear  distortion  capability  of  the  Q-19  element  are 
presented  in  Figure  8.  The  structure  is  a  circular  plate,  simply  supported  at  the  outer 
boundary  with  a  transverse  shearing  force  applied  uniformly  about  the  edge  of  a  central 
hole  (Reference  25).  Taking  advantage  of  rotational  symmetry  about  the  hole,  it  was  necessary 
to  consider  only  a  sector  of  the  plate  in  the  finite  element  idealization,  as  shown.  Comparison 
of  the  finite  element  results  with  the  exact  theory  demonstrates  the  effectiveness  of  the 
assumed  shear  distortion  mechanism. 


The  shear  stress  resultants  Q  computed  in  a  second  analysis  with  the  Q-19  shear 
distortion  element  are  shown  in  Figure  9.  The  structure  in  this  case  was  a  square,  simply 
supported  plate  subjected  to  a  uniform  loading  (Figure  9a)  and  a  central  concentrated  loading 
(Figure  9b).  The  mesh  was  8  x  8  on  the  l/4  plate  system;  the  plate  thickness  was  1/10  of  the 
total  span.  Although  no  other  solution  is  available  to  compare  with  these  results,  it  is  of 
interest  to  note  that  shear  distortion  caused  an  increase  in  central  deflection  over  the 
results  of  Figure  6  of  10%  under  uniform  loading  and  294%  with  the  concentrated  load. 
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Plate  Section 
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Finite  Element  Idealization 


Figure  8.  Analysis  of  Thick  Annular  Plate 
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Figure  9.  Computer  Plots  of  the  Transverse  Shear  Q-X 
in  a  Simply  Supported  Square  Plate 
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ANALYSIS  OF  PLATE  VIBRATIONS 
Eigenvalue  Equation 

The  equation  of  motion  of  a  discrete  coordinate  system  in  free  vibration  may  be  written: 

K  v  =  w2  M  v  (6i ) 

in  which  K  represents  the  stiffness  of  the  assembled  structure,  M  is  the  inertia  matrix 
of  the  system,  oj  is  the  circular  frequency  of  vibration  and  v  is  the  vibration  mode  shape. 
The  stiffness  matrix  used  here  is  the  same  as  would  be  used  for  a  static  analysis;  in  a 
finite  element  solution  it  is  formed  by  appropriate  superposition  of  the  individual  element 
stiffnesses.  The  mass  matrix  may  similarly  be  formed  by  assembly  of  the  individual  element 
mass  matrices.  The  first  problem  to  be  considered,  then,  is  the  definition  of  the  element 
mass  matrix. 

Element  Mass  Matrix 


Two  different  formulations  of  the  element  mass  have  been  used  extensively  in  finite 
element  studies:  the  lumped  mass  (LM)  and  the  consistent  mass  (CM)  idealizations.  For 
the  lumped  mass  approach,  it  is  assumed  that  the  mass  is  concentrated  in  points  at  the 
nodes  of  the  element,  with  nodal  values  chosen  to  be  statically  equivalent  to  the  actual 
mass  distribution.  In  the  case  of  a  uniform  triangular  element,  one  third  of  the  mass  is 
concentrated  at  each  comer.  The  LM  matrix  is  diagonal  in  form,  with  non-zero  terms 
associated  only  with  the  translational  degrees  of  freedom. 


The  CM  matrix  may  be  derived  from  the  displacement  interpolation  functions  assumed 
for  the  element.  Using  the  displacements  given  in  Equation  17,  the  consistent  mass  matrix 
for  each  subelement  is 


1 


dA 


162) 


where  p  is  the  mass  per  unit  area.  The  complete  element  CM  matrix  is  obtained  by  arming 
the  contributions  of  the  individual  subelements  (by  analogy  with  Equation  32).  The  CM  matrix 
is  fully  populated,  in  contrast  to  the  diagonal  LM  matrix. 


In  forming  the  LM  matrix  for  the  Q-19  element,  it  was  considered  desirable  to  retain 
the  mass  at  the  interior  nodal  point;  thus,  this  nodal  point  was  retained  also  in  the  stiffness 
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Degrees  of  Freedom  for  1/4  Plate: 
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Figure  10.  Reduced  Frequencies  X  mn  =  a)mn  a2^/—  of  S.S.  Square  Plate 
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Quadrilateral  Idealizations: 


Degrees  of  Freedom  for  1/2  Plate: 
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(1)  Report  DRL- 231,  Def.  Res.  Lab.,  Uni  versity  of  Texas,  Austin  ,  1949 

(2)  J.  Aero.  Sci.,  £0,  331,  1951 

(3)  J.Appl.  Mech.,18,  2,  1951 

(4)  J.ASME,  March  1962,  p.  315. 


Figure  11. 


Reduced  Frequencies  X 


of  Square  Cantilever  Plate 
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Quod ritateral  Idealizations: 


Poisson's  ratio  v  =  0.29 


Red. 
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( 1 )  Report  DRL-  231,  Def,  Res.  Lab.,  University  of  Texas,  Austin,  1949 

(2)  J.  Appl.  Mech.  ,_18,  2,  1951 


Figure  12. 


Reduced  Frequencies  X  n 
Plate  (Skewangle  =  30°) 


of  a  Rhombic  Cantilever 
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Quadrilateral  Idealizatiors: 


Degrees  of  Freedom  for  l/E  Plate: 
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Figure  13. 
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matrix  for  vibration  analysis.  The  CM  matrix  for  the  Q-19  element  was  developed  by  as¬ 
sembling  the  CM  matrices  for  4  LCCT-9  elements;  i.e.  the  interior  mid- side  nodes  were 
omitted.  The  stiffness  matrix  used  in  the  CM  vibration  studies  was  the  same  as  that  used 
in  the  LM  studies,  the  assembled  CM  matrix  for  the  complete  structure  had  a  similar  banded 
form  to  that  of  the  stiffness  matrix. 

Analytical  Results 

The  vibration  frequencies  computed  for  a  series  of  plate  systems  are  presented  in 
Figures  10-13,  The  systems  comsidered  were  a  rectangular  simply  supported  plate,  a 
rectangular  cantilever,  a  skewed  cantilever,  and  a  triangular  cantilever.  Different  mesh 
sizes  were  used  in  each  case  to  demonstrate  convergence,  and  both  LM  and  CM  idealizations 
were  utilized.  Also  shown  are  results  of  other  analyses  and  some  experimental  results.  It 
is  apparent  in  these  figures  that  the  results  of  the  finite  element  analysis  are  excellent, 
giving  reliable  agreement  with  experiments  and  with  exact  theories.  Of  greatest  significance 
is  the  comparison  between  LM  and  CM  analyses;  in  almost  every  case,  the  LM  result  for  a 
given  mesh  size  is  more  accurate  than  the  CM  value,  in  spite  of  the  fact  that  the  CM  systems 
involved  2  to  3  times  more  degrees  of  freedom.  The  only  advantage  of  the  CM  result  is  that 
it  provides  an  upper  bound  to  the  exact  value,  while  the  LM  result  may  be  either  low  or  high. 
Plots  of  a  few  mode  shapes  for  the  skewed  cantilever  plate  are  shown  in  Figure  14. 

ANALYSIS  OF  PLATE  BUCKLING 
Eigenvalue  Equation 

The  linearized  buckling  analysis  of  plates  leads  to  an  eigenvalue  equation  which  is  very 
similar  in  form  to  the  vibration  equation;  it  may  be  written: 

K  v  :  X  Kg  W  (63) 

in  which  KQ  is  the  initial  stress  (or  geometric  stiffness)  matrix  of  the  system  and  X  is 
the  critical  load  factor,  with  other  terms  as  in  Equation  59.  Since  the  elastic  stiffness 
matrix  K  used  here  is  the  same  as  that  used  for  ordinary  static  analyses,  the  essential 
problem  in  the  stability  analysis  is  evaluation  of  the  initial  stress  matrix. 

Initial  Stress  Matrix 

The  initial  stress  matrix  for  plate  buckling  analysis  may  be  interpreted  physically  as  the 
out-of-plane  nodal  forces  resulting  from  the  action  of  existing  in-plane  (membrane)  stresses 
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on  the  out-of -plane  displacements.  It  may  be  calculated,  in  general,  from  an  expression 
of  the  form 


*6  -L  [£,x  *l,y  } 


Nx  Nxy 
Nxy  Ny 


L^w.yj 


dA 


(64) 


in  which  Nx<  and  Ny  represent  the  membrane  stress  resultant  components  existing  in 
the  element,  and  <f>  w  x  y  rePresent  the  slopes  of  the  assumed  displacement 

interpolation  functions  in  the  x  and  y  directions. 


As  was  the  case  in  the  development  of  mass  matrices,  different  levels  of  approximation 
can  be  employed  in  the  formulation  of  the  initial  stress  matrices;  i.e.,  different  orders  of 
displacement  expansions  may  be  used  for  K-  The  simplest  available  approximation  is 
based  on  a  single  linear  expansion  (SL)  for  the  entire  triangle;  i.e.,  for  this  case 

<SL)  -  [  £,  C2  £,]  165) 

The  next  higher  order  of  practical  usage  is  the  single  cubic  expansion  (SC)  where  the  dis¬ 
placements  are  those  assumed  in  developing  the  incompatible  triangle  (here  denoted  BCIZ) 
of  Reference  13.  The  highest  order  considered  here  is  the  compatible  cubic  expansion,  in 
which  the  interpolation  functions  are  those  defined  for  each  subelement  in  Equation  17,  i.e.. 


(  LCCT)  * 


A(i) 

9 


(66) 


In  the  case  where  the  initial  stress  matrix  is  based  on  the  same  displacement  expansions 
used  in  deriving  the  stiffness  matrix,  it  may  be  termed  the  consistent  initial  stress  matrix; 
thus  the  consistent  K  ^  for  use  with  the  LCCT  elements  is  obtained  by  using  Equation  66 
in  Equation  64,  whereas  the  consistent  KQ  for  an  analysis  in  which  the  stiffness  matrix 
was  based  on  the  single  cubic  expansion  would  require  the  use  of  that  same  expansion  in 
Equation  64. 


Analytical  Results 


Results  of  a  series  of  buckling  analyses  of  three  different  rectangular  plates  are  presented 
in  Figure  15.  Part  “a"  shows  the  case  of  a  square  simply  supported  plate  subjected  to  uniaxial 
membrane  stress  N^;  Part  “b”  is  a  square  clamped  plate  with  uniform  biaxial  membrane 
stresses  Nx  =  Ny;  and  Part  “c”  is  a  simply  supported  rectangular  plate  (5/4  aspect  ratio) 
subjected  to  uniform  shear  stress  N  .  In  each  case,  five  different  finite  element  types 
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(  a)  S.  S.  Square  Plate  is 
Uniaxial  Compression 


4.000  12 

16 


I* —  o-*l 


(b)  Clamped  Square  Plate  in 
Bioxial  Compression 
a-  =  <r 


j  i 


w 1 


(c)  S.  S.  Rectangulor  Plate 
(h/o  =  . 8 )  under  Pure 
Shear 
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NOTES!  (*)  Only  1/4  of  the  plote  was  actually  considered. 

(*•)  These  values  (  calculated  from  energy  methods)  are  also  approximations  ,  and  the  last  digit  is 
not  guaranteed . 


Figure  15.  Buckling  Analysis  of  Three  Model  Plate  Problems 
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were  employed  for  each  of  several  different  mesh  systems.  Using  the  Q-19  stiffness,  initial 
stress  stiffnesses  were  computed  for  the  SL,  SC,  and  the  consistent  displacement  expansions. 
In  addition,  the  BCIZ  and  ACM  elements  were  used  with  their  corresponding  consistent 
initial  stress  matrices  (the  latter  results  being  from  Reference  33). 

Comparison  of  these  results  shows  that  the  Q-19  element  with  the  SC  or  consistent 
initial  stress  matrices  gave  consistently  superior  results.  The  single  linear  expansion 
K  q(SL),  although  easily  formed,  is  not  sufficiently  accurate  to  be  recommended  in  general. 
A  curious  feature  of  the  BCIZ  results  is  that  they  do  not  show  regular  improvement  with  mesh 
refinement  —  in  fact,  they  tend  to  give  poorer  results  with  the  finer  meshes  in  the  shear 
stress  loading.  In  these  analyses,  the  ACM  results  are  seen  to  be  well  behaved,  although  not 
as  accurate  as  the  Q-19  results;  in  addition,  it  must  be  remembered  that  they  also  impose 
a  significant  geometric  limitation  because  of  their  rectangular  shape. 

A  significant  conclusion  that  may  be  drawn  from  these  results  is  that  the  initial  stress 
matrix  based  on  the  single  cubic  expansion  is  as  effective  as  the  consistent  stiffness  for  the 
Q-19  element.  This  initial  stress  matrix  is  highly  recommended,  therefore,  because  it  requires 
much  less  computational  effort  to  formulate.  From  this  observation,  it  appears  that  lower 
order  expansions  may  be  used  generally  in  deriving  the  geometric  stiffness  than  are  required 
in  the  formulation  of  the  elastic  stiffness.  This  conclusion  is  parallel  to  the  observation  made 
in  the  vibration  analyses  regarding  the  relative  merits  of  the  LM  and  CM  matrices,  where 
the  lower  order  approach  also  was  considered  superior. 

The  computed  buckled  shapes  for  the  fine  mesh  Q-19- analyses  of  these  three  structural 
systems  are  shown  in  Figure  16.  Also  shown  is  the  buckled  shape  for  a  shear  loaded  rectangle 
with  an  aspect  ratio  of  5/2,  showing  the  two  lobe  deflection  pattern. 
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(A)  Simply  Supported  Square  Plate  Under  Uniform  (B)  Clamped  Square  Plate  Under  Uniform  Biaxial 
Uniaxial  Compression  Sig  -  XX  Compression 


(C)  Simply  Supported  Rectangular  Plate  (8/A  =  0.8) 
Under  Pure  Shear 


Figure  16.  Computer  Plots  of  Plate  Buckling  Modes 
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CONCLUSIONS 

The  principal  conclusions  which  may  be  drawn  from  the  study  reported  here  are  as 
follows: 

(1)  The  Q-19  fully  compatible  quadrilateral  appears  to  be  the  most  efficient  general 
plate  bending  element  yet  developed.  Because  it  has  only  12  degrees  of  freedom  when  as¬ 
sembled  into  the  complete  system,  and  no  mid-side  nodes,  it  offers  maximum  computational 
efficiency. 

(2)  The  shear  distortion  mechanism  incorporated  into  this  element  provides  a  reasonable 
approximation  of  thick  plate  behavior,  while  leaving  the  general  form  of  the  element  stiffness 
matrix  unchanged.  Thus,  it  may  easily  be  incorporated  into  existing  plate  analysis  programs. 

(3)  The  plate  vibration  studies  reported  here  demonstrate  that  a  lumped  mass  matrix  is 
more  efficient  in  representing  inertial  effects  than  a  consistent  mass  matrix.  (The  same 
relative  behavior  had  been  observed  previously  in  dynamic  axi-symmetric  analyses.) 

(4)  A  similar  conclusion  can  be  drawn  from  the  plate  buckling  results:  the  initial  stress 
matrix  based  on  a  single  cubic  expansion  is  computationally  more  efficient  than  the  higher 
order  consistent  matrix  for  the  Q-19  element. 

It  also  is  of  interest  to  report  that  the  LCCT  and  Q-19  elements  have  been  used  extensively 
and  effectively  in  providing  the  bending  stiffness  infinite  element  thin  shell  analysis  programs 
(Reference  34). 
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APPENDIX 

COMPATIBLE  DISPLACEMENT  FUNCTIONS 

For  subelement  3,  Equation  17  gives  w(3)  =  <£(3)  r  where  the  displacement  functions 
may  be  expressed:  . 

$(3)-  ^3)  ^3)£(3)  s{3)  <£131  i(3>  ^(3)^(3)  <A(3)  ^(l)^(3)At3) 

Wl  0X1  0y  I  w2  0x2  0y2  w3  0(3  0y3  04  05  06 

in  which  (in  terms  of  the  dimensions  and  coordinates  of  the  complete  element) : 

*1' ‘  C2,13-2?,'-1-  6m3C,?2C3H3  [3IX2-m5)C,+  (2^3-x2)C3-3^  c2] 


0X  I 


(b«c8-b8c2)+(b1-b^8)elclc. 


4> 


(3) 

W2 


+  6  [3<b*  X2  +  b3  /t3-2l>ll£,  +3'*>3  >CZ  +  t3b,  -t>2X2-2bj/is)C3] 

*  C2(3-2C2'+  6X3c,C2C3+C3  [3<f,-x3)C2+12V/VC3-3X3c,  ] 


$<•> 

0x2 


c;(63Crb,c2H-(t3  x3-b2t  c,c2c3 


+  ?^[3(2b2  -b3X3-b,ft»  C2+3lb2-b3X3)C1+l-3b2-blM,«b3XJl  £.] 

=  t3[3l'^2^,<-3"  +  x,  >C,+  m2-x,)C3] 

=  ■g'  [3(3bi+b2,_b.Xi>C2  +  ‘b>z  /“2-b,  X,)C3-  3  (b,+3b2  +b2  /i2)C,  ] 

■■S£K t.i,<  <5c3-3'] 

•^-[c:<3c2c3i] 


<£<*> 

‘  W3 

(9x3 

a!51 

^04 
A  (3) 

^05 

$r  ■  $  towi 


3L, 


For  (i  =  1,2,3)  change  all  b’s  in  <^3*  to  a’s. 

0yi  0xi 

Above  expressions  apply  to  subelement  3,  where  £(  >  £^  ,  £2  >  £ 

In  subelements  1  (where  £  >  £  ,  £  >£  )  and  in  subelement  2  (where  £  >  £  £  >£  )  permute 

2  13  1  3  2  12 

cyclically  all  subscripts  and  superscripts  (1-2-3  permutes  to  2-3-1  and  4-5-6  to  5-6-4); 
for  example:  =  t(  2  [  3  (l  +/x3 )  £2  +  3  (  I  +  )  £3  +  ( |  -  -  \gj  £(  ]  . 


440 


